execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

ex8

non linear regression 
with one explanatory variable  
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}
single linear regression

ex8-0.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)
## Initial log joint probability = -594.909 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24      -39.8541   0.000404967   0.000456077           1           1       34    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -39.85
##    b0         4.87
##    b1         1.84
##    s          4.45
##    m0[1]     11.26
##    m0[2]     18.67
##    m0[3]      9.94
##    m0[4]     24.60
##    m0[5]     32.34
##    m0[6]     24.19
##    m0[7]     31.27
##    m0[8]     11.26
##    m0[9]      9.06
##    m0[10]    31.84
##    m0[11]    27.11
##    m0[12]    30.40
##    m0[13]     9.48
##    m0[14]    35.05
##    m0[15]    36.03
##    m0[16]     4.93
##    m0[17]    28.48
##    m0[18]    33.52
##    m0[19]    15.69
##    m0[20]     5.45
##    y0[1]     15.84
##    y0[2]     12.43
##    y0[3]      8.42
##    y0[4]     22.29
##    y0[5]     34.28
##    y0[6]     24.69
##    y0[7]     31.86
##    y0[8]     15.69
##    y0[9]      0.22
##    y0[10]    30.34
##    y0[11]    31.39
##    y0[12]    29.04
##    y0[13]     2.06
##    y0[14]    35.01
##    y0[15]    30.44
##    y0[16]     6.25
##    y0[17]    28.82
##    y0[18]    35.81
##    y0[19]    19.67
##    y0[20]     2.41
##    m1[1]      4.93
##    m1[2]      8.39
##    m1[3]     11.84
##    m1[4]     15.30
##    m1[5]     18.75
##    m1[6]     22.21
##    m1[7]     25.66
##    m1[8]     29.12
##    m1[9]     32.57
##    m1[10]    36.03
##    y1[1]      3.82
##    y1[2]     10.88
##    y1[3]     13.01
##    y1[4]      5.94
##    y1[5]     13.24
##    y1[6]     18.10
##    y1[7]     27.49
##    y1[8]     25.23
##    y1[9]     36.65
##    y1[10]    31.63
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -39.95 -39.65 1.23 1.07 -42.36 -38.53 1.00      535      991
##    b0       4.96   4.89 2.05 2.15   1.50   8.32 1.00      251      276
##    b1       1.83   1.84 0.19 0.19   1.52   2.15 1.00      369      724
##    s        5.07   4.93 0.92 0.84   3.78   6.73 1.00     1161     1242
##    m0[1]   11.31  11.27 1.55 1.62   8.73  13.90 1.00      265      309
##    m0[2]   18.68  18.66 1.19 1.18  16.75  20.67 1.00      429      816
##    m0[3]    9.99   9.94 1.64 1.73   7.27  12.71 1.00      259      262
##    m0[4]   24.57  24.57 1.21 1.16  22.57  26.52 1.00     1779     1353
##    m0[5]   32.27  32.23 1.63 1.56  29.62  35.01 1.00     2559     1522
##    m0[6]   24.16  24.16 1.20 1.15  22.21  26.10 1.00     1574     1358
##    m0[7]   31.21  31.17 1.56 1.48  28.68  33.81 1.00     2641     1430
##    m0[8]   11.31  11.27 1.55 1.62   8.73  13.90 1.00      265      309
##    m0[9]    9.12   9.05 1.71 1.80   6.29  11.93 1.00      256      250
##    m0[10]  31.77  31.73 1.60 1.53  29.17  34.47 1.00     2597     1498
##    m0[11]  27.07  27.05 1.31 1.23  24.95  29.21 1.00     2576     1386
##    m0[12]  30.35  30.32 1.50 1.43  27.92  32.83 1.00     2684     1368
##    m0[13]   9.54   9.48 1.68 1.77   6.76  12.31 1.00      257      258
##    m0[14]  34.96  34.96 1.84 1.75  31.96  38.06 1.00     2204     1471
##    m0[15]  35.94  35.94 1.92 1.84  32.76  39.14 1.00     2050     1366
##    m0[16]   5.01   4.95 2.04 2.15   1.56   8.36 1.00      251      276
##    m0[17]  28.44  28.41 1.38 1.29  26.18  30.71 1.00     2689     1409
##    m0[18]  33.44  33.42 1.72 1.65  30.66  36.31 1.00     2442     1451
##    m0[19]  15.71  15.68 1.29 1.28  13.61  17.85 1.00      320      588
##    m0[20]   5.53   5.47 2.00 2.11   2.16   8.80 1.00      251      267
##    y0[1]   11.17  11.03 5.38 4.81   2.19  20.42 1.00     1395     1595
##    y0[2]   18.65  18.72 5.28 5.17  10.29  26.93 1.00     1921     1891
##    y0[3]    9.96   9.99 5.43 5.33   0.94  18.67 1.00     1322     1839
##    y0[4]   24.53  24.53 5.40 5.18  15.55  33.39 1.00     1927     1821
##    y0[5]   32.10  32.08 5.42 5.08  23.16  40.95 1.00     1950     1838
##    y0[6]   24.16  24.16 5.42 5.14  15.17  32.97 1.00     1936     1787
##    y0[7]   31.35  31.42 5.33 5.05  22.65  40.07 1.00     2099     1932
##    y0[8]   11.45  11.47 5.37 5.05   2.76  20.16 1.00     1257     1731
##    y0[9]    8.97   8.97 5.49 5.18   0.08  18.11 1.00     1448     1573
##    y0[10]  31.63  31.69 5.51 5.35  22.71  40.56 1.00     2121     1856
##    y0[11]  27.19  27.24 5.42 5.14  18.31  35.90 1.00     1953     1999
##    y0[12]  30.31  30.29 5.50 5.18  21.23  39.30 1.00     1955     1913
##    y0[13]   9.47   9.41 5.39 5.09   0.68  18.47 1.00     1282     1598
##    y0[14]  34.95  34.92 5.29 5.10  26.67  43.65 1.00     2148     1685
##    y0[15]  36.05  36.00 5.60 5.21  26.94  45.43 1.00     2045     1855
##    y0[16]   5.05   4.86 5.49 5.20  -3.86  14.48 1.00      964     1464
##    y0[17]  28.41  28.46 5.31 5.27  19.51  37.15 1.00     2130     1799
##    y0[18]  33.22  33.23 5.35 5.13  24.73  41.92 1.00     2271     1889
##    y0[19]  15.58  15.60 5.35 5.24   6.78  24.24 1.00     1752     1836
##    y0[20]   5.49   5.43 5.62 5.59  -3.84  14.55 1.01     1187     1796
##    m1[1]    5.01   4.95 2.04 2.15   1.56   8.36 1.00      251      276
##    m1[2]    8.45   8.39 1.76 1.84   5.55  11.37 1.00      254      278
##    m1[3]   11.89  11.84 1.51 1.57   9.37  14.43 1.00      268      345
##    m1[4]   15.32  15.28 1.31 1.31  13.18  17.49 1.00      312      531
##    m1[5]   18.76  18.75 1.19 1.18  16.84  20.74 1.00      434      816
##    m1[6]   22.19  22.18 1.16 1.14  20.32  24.11 1.00      835     1059
##    m1[7]   25.63  25.62 1.25 1.17  23.59  27.65 1.00     2174     1357
##    m1[8]   29.06  29.04 1.42 1.33  26.77  31.39 1.00     2706     1409
##    m1[9]   32.50  32.46 1.65 1.58  29.83  35.27 1.00     2539     1523
##    m1[10]  35.94  35.94 1.92 1.84  32.76  39.14 1.00     2050     1366
##    y1[1]    4.96   5.00 5.52 5.28  -4.41  14.03 1.00     1097     1551
##    y1[2]    8.36   8.46 5.32 5.11  -0.51  17.23 1.00     1269     1732
##    y1[3]   12.04  12.05 5.11 4.82   4.06  20.46 1.00     1329     1867
##    y1[4]   15.27  15.28 5.06 4.84   7.03  23.52 1.00     1971     1894
##    y1[5]   18.73  18.83 5.47 5.20   9.61  27.48 1.00     1954     1934
##    y1[6]   22.39  22.41 5.26 5.04  13.87  30.75 1.00     1603     1890
##    y1[7]   25.59  25.71 5.39 5.29  16.60  34.34 1.00     1950     1891
##    y1[8]   29.18  29.33 5.22 4.93  20.54  37.91 1.00     2190     1992
##    y1[9]   32.62  32.80 5.31 5.19  23.95  40.86 1.00     2106     1499
##    y1[10]  35.98  35.96 5.52 5.26  27.06  45.20 1.00     2011     1941

quadratic regression  
y=b0+b2(x-b1)**2  

ex8-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real b2;
  real<lower=0> s;
}
model{
  y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b2*(x[i]-b1)^2;
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b2*(x1[i]-b1)^2;
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)
## Initial log joint probability = -316802 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -35.4836       6.06152       279.657           1           1      199    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199      -35.4468      0.686777       1461.74      0.2051      0.2051      348    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      299      -35.4432       5.73746       1968.49           1           1      495    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      399       -35.442       0.46783       2597.15      0.3871      0.3871      635    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      499      -35.4412      0.666778       5867.63      0.3691      0.3691      782    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      599      -35.4408        4.5187       4380.86           1           1      934    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      699      -35.4405     0.0200503       262.313       0.153      0.8623     1082    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      799      -35.4403       5.26346       9198.39       0.566           1     1235    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      831      -35.4403       1.52194       625.731           1           1     1276    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -35.44
##    b0       166.07
##    b1      2585.60
##    b2         0.00
##    s          3.56
##    m0[1]      9.09
##    m0[2]      8.67
##    m0[3]      8.72
##    m0[4]      8.13
##    m0[5]      8.80
##    m0[6]      8.75
##    m0[7]      8.14
##    m0[8]      9.06
##    m0[9]      8.39
##    m0[10]     8.41
##    m0[11]     9.07
##    m0[12]     9.16
##    m0[13]     8.93
##    m0[14]     8.22
##    m0[15]     8.20
##    m0[16]     8.88
##    m0[17]     8.82
##    m0[18]     8.33
##    m0[19]     8.26
##    m0[20]     8.47
##    y0[1]     10.27
##    y0[2]     10.05
##    y0[3]      9.86
##    y0[4]      0.05
##    y0[5]      6.96
##    y0[6]     12.64
##    y0[7]     10.14
##    y0[8]      3.26
##    y0[9]      5.32
##    y0[10]     6.55
##    y0[11]     7.04
##    y0[12]     7.29
##    y0[13]     7.58
##    y0[14]     5.47
##    y0[15]     8.17
##    y0[16]     9.53
##    y0[17]    15.14
##    y0[18]     6.79
##    y0[19]    11.42
##    y0[20]    11.09
##    m1[1]      8.13
##    m1[2]      8.25
##    m1[3]      8.36
##    m1[4]      8.47
##    m1[5]      8.59
##    m1[6]      8.70
##    m1[7]      8.82
##    m1[8]      8.93
##    m1[9]      9.04
##    m1[10]     9.16
##    y1[1]     15.09
##    y1[2]      7.19
##    y1[3]      7.32
##    y1[4]     12.51
##    y1[5]      6.45
##    y1[6]     11.89
##    y1[7]      9.87
##    y1[8]      6.61
##    y1[9]      7.75
##    y1[10]     8.60
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.31 -11.97 1.49 1.23 -15.25 -10.59 1.01      530      461
##    b0       4.21   4.20 0.41 0.39   3.55   4.88 1.01      384      498
##    b1       3.96   3.96 0.09 0.08   3.82   4.11 1.00     1931      992
##    b2       0.56   0.56 0.04 0.04   0.48   0.63 1.01      548      742
##    s        1.18   1.15 0.21 0.19   0.88   1.57 1.00      917      866
##    m0[1]   13.87  13.90 0.56 0.54  12.92  14.76 1.00     1897     1414
##    m0[2]    4.45   4.45 0.41 0.39   3.80   5.12 1.01      391      475
##    m0[3]    4.85   4.84 0.40 0.37   4.21   5.50 1.01      407      554
##    m0[4]   11.92  11.91 0.59 0.57  10.97  12.86 1.00     1912     1464
##    m0[5]    5.99   5.99 0.37 0.33   5.40   6.60 1.01      477      718
##    m0[6]    5.26   5.25 0.39 0.35   4.64   5.89 1.01      427      601
##    m0[7]   11.70  11.69 0.57 0.56  10.78  12.62 1.00     1952     1513
##    m0[8]   12.78  12.80 0.50 0.49  11.94  13.57 1.00     2069     1445
##    m0[9]    5.70   5.69 0.35 0.33   5.12   6.28 1.01      514      686
##    m0[10]   5.36   5.35 0.36 0.34   4.77   5.95 1.01      466      643
##    m0[11]  13.08  13.11 0.52 0.51  12.22  13.91 1.00     2037     1458
##    m0[12]  16.51  16.52 0.73 0.70  15.27  17.66 1.00     1276     1250
##    m0[13]   8.66   8.67 0.35 0.35   8.09   9.23 1.00     1055      882
##    m0[14]   9.28   9.28 0.43 0.41   8.61   9.97 1.00     1949     1403
##    m0[15]   9.84   9.83 0.46 0.44   9.11  10.59 1.00     2000     1567
##    m0[16]   7.36   7.36 0.35 0.32   6.79   7.92 1.01      656      802
##    m0[17]   6.27   6.27 0.36 0.33   5.68   6.86 1.01      503      773
##    m0[18]   6.62   6.62 0.35 0.33   6.05   7.18 1.00      728      820
##    m0[19]   8.28   8.28 0.38 0.36   7.67   8.91 1.00     1591     1060
##    m0[20]   4.68   4.67 0.38 0.36   4.05   5.32 1.01      406      500
##    y0[1]   13.87  13.89 1.33 1.26  11.72  16.06 1.00     2109     1654
##    y0[2]    4.50   4.48 1.28 1.19   2.38   6.65 1.00     1591     1632
##    y0[3]    4.81   4.84 1.28 1.22   2.76   6.91 1.00     1680     1775
##    y0[4]   11.89  11.91 1.34 1.30   9.72  14.09 1.00     2011     1711
##    y0[5]    5.95   5.94 1.27 1.22   3.85   8.03 1.00     1787     1886
##    y0[6]    5.25   5.25 1.24 1.26   3.14   7.29 1.00     1489     1312
##    y0[7]   11.75  11.75 1.36 1.30   9.53  13.90 1.00     1546     1593
##    y0[8]   12.72  12.71 1.34 1.25  10.48  14.90 1.00     2006     1858
##    y0[9]    5.68   5.67 1.29 1.23   3.58   7.78 1.00     1739     1536
##    y0[10]   5.37   5.35 1.24 1.22   3.29   7.40 1.00     1874     1819
##    y0[11]  13.11  13.10 1.31 1.22  10.97  15.32 1.00     1863     1605
##    y0[12]  16.48  16.51 1.36 1.32  14.22  18.61 1.00     1674     1631
##    y0[13]   8.66   8.67 1.23 1.20   6.70  10.74 1.00     2028     1878
##    y0[14]   9.27   9.25 1.25 1.16   7.21  11.39 1.00     1901     1656
##    y0[15]   9.81   9.80 1.27 1.23   7.79  11.89 1.00     2030     1932
##    y0[16]   7.36   7.35 1.27 1.23   5.27   9.46 1.00     1861     1860
##    y0[17]   6.30   6.30 1.25 1.17   4.26   8.36 1.00     1632     1858
##    y0[18]   6.61   6.63 1.24 1.15   4.49   8.61 1.00     1834     1767
##    y0[19]   8.31   8.29 1.26 1.22   6.24  10.41 1.00     1815     1839
##    y0[20]   4.67   4.63 1.25 1.17   2.68   6.73 1.00     1673     1765
##    m1[1]   11.92  11.91 0.59 0.57  10.97  12.86 1.00     1912     1464
##    m1[2]    8.53   8.53 0.39 0.38   7.90   9.18 1.00     1710     1179
##    m1[3]    6.12   6.11 0.35 0.33   5.55   6.69 1.00      593      759
##    m1[4]    4.68   4.67 0.38 0.36   4.05   5.32 1.01      406      500
##    m1[5]    4.22   4.21 0.41 0.39   3.56   4.89 1.01      383      501
##    m1[6]    4.73   4.72 0.40 0.38   4.08   5.38 1.01      402      528
##    m1[7]    6.21   6.21 0.36 0.33   5.62   6.80 1.01      497      768
##    m1[8]    8.67   8.67 0.35 0.35   8.09   9.23 1.00     1057      882
##    m1[9]   12.10  12.12 0.47 0.46  11.30  12.83 1.00     2115     1425
##    m1[10]  16.51  16.52 0.73 0.70  15.27  17.66 1.00     1276     1250
##    y1[1]   11.92  11.91 1.36 1.29   9.79  14.18 1.00     2051     1969
##    y1[2]    8.55   8.57 1.28 1.21   6.48  10.67 1.00     1875     1854
##    y1[3]    6.16   6.18 1.25 1.18   4.06   8.23 1.00     1529     1698
##    y1[4]    4.65   4.63 1.21 1.14   2.69   6.63 1.00     1304     1513
##    y1[5]    4.15   4.16 1.24 1.23   2.10   6.22 1.00     1417     1736
##    y1[6]    4.66   4.68 1.27 1.22   2.56   6.72 1.00     1351     1891
##    y1[7]    6.20   6.19 1.28 1.21   4.04   8.27 1.00     1646     1682
##    y1[8]    8.65   8.66 1.28 1.29   6.56  10.72 1.00     1724     1728
##    y1[9]   12.07  12.10 1.30 1.23   9.94  14.18 1.00     1916     1677
##    y1[10]  16.51  16.55 1.41 1.36  14.12  18.77 1.00     1530     1722

both log regression  
log y=b0+b1*log x   x,y>0
y=exp b0 * x**b1

ex8-2.stan

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector<lower=0>[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~lognormal(b0+b1*log(x),s);
}
generated quantities{
  vector[N] m0;
  vector<lower=0>[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*log(x[i]);
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector<lower=0>[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*log(x1[i]);
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
             qplot(log(x),log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -67.6061 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -30.3625   0.000408541   0.000581496           1           1       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -30.36
##    b0         0.78
##    b1         2.22
##    s          1.10
##    m0[1]      5.64
##    m0[2]      0.80
##    m0[3]      5.07
##    m0[4]      4.71
##    m0[5]      4.77
##    m0[6]      4.63
##    m0[7]      5.65
##    m0[8]      4.10
##    m0[9]      4.29
##    m0[10]    -0.21
##    m0[11]     4.80
##    m0[12]     4.77
##    m0[13]     5.57
##    m0[14]     4.81
##    m0[15]     4.52
##    m0[16]     3.86
##    m0[17]     5.63
##    m0[18]     3.20
##    m0[19]     2.06
##    m0[20]    -0.75
##    y0[1]     78.58
##    y0[2]      0.59
##    y0[3]     46.42
##    y0[4]    301.26
##    y0[5]    122.24
##    y0[6]     25.20
##    y0[7]    523.62
##    y0[8]    114.82
##    y0[9]    229.40
##    y0[10]     2.21
##    y0[11]    61.53
##    y0[12]    29.10
##    y0[13]   318.82
##    y0[14]   161.48
##    y0[15]    13.98
##    y0[16]   152.85
##    y0[17]   164.69
##    y0[18]    36.11
##    y0[19]     6.96
##    y0[20]     0.82
##    m1[1]     -0.75
##    m1[2]      1.59
##    m1[3]      2.70
##    m1[4]      3.44
##    m1[5]      4.00
##    m1[6]      4.44
##    m1[7]      4.81
##    m1[8]      5.13
##    m1[9]      5.40
##    m1[10]     5.65
##    y1[1]      0.56
##    y1[2]     40.36
##    y1[3]     30.99
##    y1[4]      9.95
##    y1[5]     80.28
##    y1[6]      7.41
##    y1[7]    150.57
##    y1[8]     48.11
##    y1[9]    224.98
##    y1[10]   180.53
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median      sd    mad     q5     q95 rhat ess_bulk ess_tail
##    lp__   -31.89 -31.54    1.33   1.05 -34.56  -30.47 1.01      772      934
##    b0       0.81   0.80    0.56   0.52  -0.08    1.75 1.01      466      545
##    b1       2.21   2.21    0.34   0.32   1.63    2.76 1.01      522      729
##    s        1.26   1.23    0.23   0.22   0.94    1.68 1.00     1635     1364
##    m0[1]    5.64   5.64    0.39   0.37   5.03    6.26 1.00     1684     1289
##    m0[2]    0.83   0.82    0.56   0.52  -0.06    1.77 1.01      466      545
##    m0[3]    5.07   5.06    0.34   0.32   4.52    5.62 1.00     2109     1384
##    m0[4]    4.71   4.71    0.31   0.29   4.21    5.22 1.00     2252     1409
##    m0[5]    4.77   4.77    0.31   0.29   4.26    5.29 1.00     2243     1330
##    m0[6]    4.63   4.63    0.31   0.28   4.13    5.13 1.00     2253     1430
##    m0[7]    5.64   5.64    0.39   0.37   5.03    6.26 1.00     1680     1289
##    m0[8]    4.11   4.11    0.29   0.27   3.64    4.57 1.00     1865     1274
##    m0[9]    4.29   4.29    0.29   0.27   3.82    4.76 1.00     2110     1381
##    m0[10]  -0.17  -0.17    0.70   0.66  -1.28    1.03 1.01      458      544
##    m0[11]   4.80   4.80    0.32   0.30   4.28    5.33 1.00     2233     1330
##    m0[12]   4.77   4.77    0.31   0.29   4.26    5.30 1.00     2242     1330
##    m0[13]   5.56   5.56    0.38   0.36   4.96    6.17 1.00     1745     1307
##    m0[14]   4.81   4.81    0.32   0.30   4.29    5.34 1.00     2232     1327
##    m0[15]   4.52   4.52    0.30   0.28   4.04    5.01 1.00     2236     1409
##    m0[16]   3.87   3.87    0.28   0.27   3.40    4.32 1.00     1451     1307
##    m0[17]   5.63   5.63    0.39   0.37   5.02    6.25 1.00     1692     1271
##    m0[18]   3.21   3.21    0.31   0.29   2.72    3.70 1.01      808     1106
##    m0[19]   2.08   2.07    0.40   0.39   1.45    2.77 1.01      520      613
##    m0[20]  -0.71  -0.72    0.78   0.74  -1.96    0.63 1.01      457      584
##    y0[1]  648.83 287.04 1401.59 316.61  28.13 2264.01 1.00     1918     1813
##    y0[2]    8.38   2.37   77.23   2.57   0.22   23.92 1.00     1335     1809
##    y0[3]  411.92 158.35 1080.41 167.01  18.81 1366.86 1.00     2132     1800
##    y0[4]  258.56 111.35  559.60 117.38  12.38  915.14 1.00     2047     1730
##    y0[5]  305.95 116.21 1127.08 118.67  14.15  936.63 1.00     1731     1878
##    y0[6]  262.61  95.38  776.56 100.43  10.02  817.51 1.00     1909     1806
##    y0[7]  743.30 282.54 2518.92 290.81  33.42 2736.99 1.00     1962     1852
##    y0[8]  144.97  64.64  280.36  67.75   7.29  530.46 1.00     2006     1837
##    y0[9]  198.00  77.98  508.93  83.65   8.62  703.79 1.00     1858     1786
##    y0[10]   2.46   0.79    6.74   0.91   0.07    8.91 1.01     1133     1666
##    y0[11] 295.55 118.52  833.31 123.40  14.02  965.03 1.00     2074     1816
##    y0[12] 344.98 118.19 1333.30 126.62  13.30 1140.70 1.00     1899     1931
##    y0[13] 664.06 265.04 1702.03 278.38  28.06 2344.33 1.00     2056     1769
##    y0[14] 342.80 127.27 1005.87 134.11  14.15 1068.39 1.00     2010     1655
##    y0[15] 251.26  86.80  984.18  90.92   9.76  871.85 1.00     2061     1684
##    y0[16] 117.42  48.66  255.17  50.69   5.80  415.75 1.00     2056     2014
##    y0[17] 796.30 274.63 3345.23 297.55  32.38 2686.79 1.00     1704     1641
##    y0[18]  62.12  24.39  203.79  25.01   2.50  190.41 1.00     2003     1661
##    y0[19]  20.85   7.62   48.99   7.96   0.85   77.17 1.00     1850     1931
##    y0[20]   1.78   0.48    6.88   0.53   0.04    6.12 1.00     1184     1676
##    m1[1]   -0.71  -0.72    0.78   0.74  -1.96    0.63 1.01      457      584
##    m1[2]    1.61   1.61    0.46   0.43   0.89    2.39 1.01      488      551
##    m1[3]    2.72   2.71    0.34   0.33   2.18    3.29 1.01      614      846
##    m1[4]    3.45   3.45    0.29   0.28   2.98    3.92 1.00      971     1211
##    m1[5]    4.00   4.01    0.29   0.27   3.54    4.46 1.00     1678     1273
##    m1[6]    4.44   4.45    0.30   0.27   3.97    4.92 1.00     2203     1388
##    m1[7]    4.81   4.81    0.32   0.30   4.29    5.34 1.00     2231     1327
##    m1[8]    5.12   5.12    0.34   0.32   4.57    5.68 1.00     2185     1384
##    m1[9]    5.40   5.40    0.37   0.34   4.82    5.98 1.00     1903     1287
##    m1[10]   5.64   5.64    0.39   0.37   5.03    6.26 1.00     1680     1289
##    y1[1]    1.68   0.48    6.20   0.55   0.04    6.36 1.00     1188     1599
##    y1[2]   13.45   4.77   69.66   5.28   0.55   41.58 1.00     1570     1724
##    y1[3]   40.28  14.95  187.71  15.60   1.72  123.67 1.00     1918     1965
##    y1[4]  103.56  32.31 1072.77  34.91   3.54  286.96 1.00     1898     1845
##    y1[5]  146.30  54.72  399.81  55.32   7.32  465.89 1.00     2187     1694
##    y1[6]  208.31  84.58  486.90  90.75   8.83  737.84 1.00     2036     1947
##    y1[7]  313.42 124.30  877.99 132.47  12.99 1136.51 1.00     1899     1855
##    y1[8]  421.69 167.36 1030.78 182.31  16.79 1459.20 1.00     2113     1968
##    y1[9]  563.86 211.92 1852.93 224.37  24.62 1864.59 1.00     1990     1877
##    y1[10] 778.48 279.18 2238.00 299.70  31.75 2664.62 1.00     1892     1766

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

qplot(log(x),log(y),col=I('red'))+
  geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=m),xy,col='black')

ex8-3

exponential increasing/decreasing  

y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)  
log y=log b0+b1* x  -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)  
b1>0 x->Infinity,y->Infinity  
b1<0 x->Infinity,y->+0  
n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)
y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)  

ex8-3-1.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0*exp(b1*x),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*exp(b1*x[i]);
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*exp(b1*x1[i]);
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)
## Initial log joint probability = -94.0922 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19         9.308   7.50058e-05   0.000441256           1           1       40    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__       9.31
##    b0         8.95
##    b1        -2.06
##    s          0.38
##    m0[1]      0.01
##    m0[2]      0.01
##    m0[3]      0.01
##    m0[4]      0.01
##    m0[5]      0.00
##    m0[6]      3.26
##    m0[7]      1.63
##    m0[8]      0.04
##    m0[9]      0.00
##    m0[10]     0.00
##    m0[11]     1.20
##    m0[12]     3.93
##    m0[13]     0.00
##    m0[14]     0.01
##    m0[15]     2.13
##    m0[16]     1.68
##    m0[17]     0.00
##    m0[18]     1.43
##    m0[19]     0.39
##    m0[20]     0.00
##    y0[1]      0.04
##    y0[2]     -0.83
##    y0[3]      0.01
##    y0[4]      0.63
##    y0[5]      0.53
##    y0[6]      3.40
##    y0[7]      1.38
##    y0[8]     -0.68
##    y0[9]      0.29
##    y0[10]    -0.45
##    y0[11]     1.52
##    y0[12]     4.11
##    y0[13]     0.03
##    y0[14]    -0.69
##    y0[15]     1.56
##    y0[16]     1.63
##    y0[17]    -0.35
##    y0[18]     1.77
##    y0[19]     0.34
##    y0[20]     0.15
##    m1[1]      3.93
##    m1[2]      1.46
##    m1[3]      0.54
##    m1[4]      0.20
##    m1[5]      0.07
##    m1[6]      0.03
##    m1[7]      0.01
##    m1[8]      0.00
##    m1[9]      0.00
##    m1[10]     0.00
##    y1[1]      4.18
##    y1[2]      1.58
##    y1[3]      0.19
##    y1[4]      0.05
##    y1[5]      0.79
##    y1[6]     -0.32
##    y1[7]      0.29
##    y1[8]     -0.39
##    y1[9]      0.03
##    y1[10]    -0.40
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__   10.51  10.83 1.31 1.14  7.94 11.95 1.01      634      781
##    b0      9.66   9.41 2.14 1.65  6.93 13.44 1.00      704      716
##    b1     -2.18  -2.14 0.35 0.31 -2.79 -1.66 1.00      751      874
##    s       0.43   0.42 0.08 0.08  0.32  0.58 1.01      801      867
##    m0[1]   0.01   0.01 0.01 0.01  0.00  0.03 1.00      780      870
##    m0[2]   0.01   0.01 0.01 0.00  0.00  0.02 1.00      778      885
##    m0[3]   0.02   0.01 0.01 0.01  0.00  0.04 1.00      781      870
##    m0[4]   0.01   0.00 0.01 0.00  0.00  0.02 1.00      778      885
##    m0[5]   0.00   0.00 0.00 0.00  0.00  0.01 1.00      772      844
##    m0[6]   3.27   3.26 0.24 0.23  2.87  3.67 1.00     1363     1181
##    m0[7]   1.58   1.58 0.19 0.18  1.27  1.88 1.00     1215     1123
##    m0[8]   0.04   0.03 0.03 0.02  0.01  0.09 1.00      790      891
##    m0[9]   0.00   0.00 0.00 0.00  0.00  0.00 1.00      770      843
##    m0[10]  0.00   0.00 0.00 0.00  0.00  0.01 1.00      774      885
##    m0[11]  1.15   1.14 0.18 0.17  0.85  1.44 1.00     1038      998
##    m0[12]  3.98   3.96 0.35 0.32  3.44  4.58 1.00      986     1137
##    m0[13]  0.00   0.00 0.00 0.00  0.00  0.00 1.00      771      843
##    m0[14]  0.01   0.01 0.01 0.01  0.00  0.04 1.00      781      870
##    m0[15]  2.09   2.08 0.19 0.17  1.78  2.39 1.00     1546     1112
##    m0[16]  1.63   1.63 0.19 0.18  1.32  1.93 1.00     1243     1122
##    m0[17]  0.00   0.00 0.00 0.00  0.00  0.00 1.00      771      843
##    m0[18]  1.38   1.38 0.19 0.18  1.07  1.68 1.00     1122     1058
##    m0[19]  0.37   0.36 0.12 0.11  0.19  0.57 1.00      865      815
##    m0[20]  0.00   0.00 0.00 0.00  0.00  0.01 1.00      774      885
##    y0[1]   0.02   0.03 0.45 0.44 -0.71  0.75 1.00     1832     1777
##    y0[2]   0.01   0.01 0.44 0.42 -0.71  0.73 1.00     1928     1799
##    y0[3]   0.03   0.04 0.44 0.43 -0.69  0.74 1.00     1946     1675
##    y0[4]   0.01   0.01 0.44 0.43 -0.71  0.76 1.00     1685     1831
##    y0[5]  -0.02  -0.02 0.44 0.41 -0.75  0.70 1.00     1960     1776
##    y0[6]   3.28   3.28 0.51 0.49  2.47  4.10 1.00     1646     1506
##    y0[7]   1.59   1.60 0.48 0.46  0.79  2.36 1.00     1935     1847
##    y0[8]   0.05   0.06 0.44 0.44 -0.67  0.74 1.00     1914     1673
##    y0[9]   0.01   0.01 0.45 0.42 -0.70  0.77 1.00     2005     1872
##    y0[10]  0.00   0.01 0.44 0.42 -0.72  0.71 1.00     1866     1703
##    y0[11]  1.15   1.15 0.49 0.46  0.35  1.94 1.00     1815     1771
##    y0[12]  3.96   3.96 0.56 0.55  3.10  4.89 1.00     1399     1178
##    y0[13] -0.01  -0.01 0.43 0.40 -0.71  0.68 1.00     1791     1870
##    y0[14]  0.02   0.03 0.43 0.40 -0.68  0.73 1.00     1920     1768
##    y0[15]  2.07   2.09 0.47 0.43  1.26  2.84 1.00     1927     1730
##    y0[16]  1.61   1.61 0.48 0.44  0.83  2.40 1.00     1800     1846
##    y0[17]  0.00   0.00 0.43 0.42 -0.73  0.69 1.00     2048     1894
##    y0[18]  1.38   1.37 0.47 0.43  0.62  2.18 1.00     1924     1723
##    y0[19]  0.36   0.36 0.46 0.43 -0.38  1.11 1.00     1835     1684
##    y0[20] -0.01  -0.01 0.44 0.42 -0.74  0.69 1.00     2024     1814
##    m1[1]   3.98   3.96 0.35 0.32  3.44  4.58 1.00      986     1137
##    m1[2]   1.40   1.40 0.19 0.18  1.10  1.70 1.00     1133     1058
##    m1[3]   0.51   0.50 0.14 0.13  0.29  0.75 1.00      889      882
##    m1[4]   0.19   0.18 0.08 0.07  0.08  0.33 1.00      824      811
##    m1[5]   0.07   0.06 0.04 0.04  0.02  0.15 1.00      798      891
##    m1[6]   0.03   0.02 0.02 0.02  0.01  0.07 1.00      787      902
##    m1[7]   0.01   0.01 0.01 0.01  0.00  0.03 1.00      780      870
##    m1[8]   0.00   0.00 0.01 0.00  0.00  0.01 1.00      775      885
##    m1[9]   0.00   0.00 0.00 0.00  0.00  0.01 1.00      772      885
##    m1[10]  0.00   0.00 0.00 0.00  0.00  0.00 1.00      771      843
##    y1[1]   3.98   3.96 0.57 0.54  3.07  4.93 1.00     1374     1467
##    y1[2]   1.41   1.40 0.48 0.46  0.61  2.22 1.00     1661     1645
##    y1[3]   0.50   0.50 0.46 0.44 -0.26  1.25 1.00     1784     1861
##    y1[4]   0.19   0.18 0.45 0.42 -0.53  0.92 1.00     1954     1755
##    y1[5]   0.08   0.08 0.44 0.43 -0.64  0.83 1.00     2001     1771
##    y1[6]   0.03   0.03 0.44 0.41 -0.68  0.76 1.00     2223     1849
##    y1[7]   0.00  -0.01 0.43 0.41 -0.69  0.71 1.00     1926     1591
##    y1[8]   0.00   0.01 0.42 0.40 -0.71  0.68 1.00     2282     1942
##    y1[9]  -0.02  -0.02 0.43 0.41 -0.75  0.70 1.00     2001     1638
##    y1[10]  0.00   0.01 0.44 0.43 -0.75  0.74 1.00     1928     1838

log y=log b0+b1* x  -> y~logN(log b0 +b1*x,s)

ex8-3-2.stan

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real<lower=-10,upper=10> b1;
  real<lower=0> s;
}
model{
  y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=log(b0)+b1*x[i];
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=log(b0)+b1*x1[i];
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -17976.8 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       48      -9.81683   1.26169e-05   0.000321397      0.2257      0.8647      135    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -9.82
##    b0         9.53
##    b1        -1.97
##    s          0.40
##    m0[1]      1.72
##    m0[2]     -1.69
##    m0[3]     -4.29
##    m0[4]      2.02
##    m0[5]      1.89
##    m0[6]     -2.69
##    m0[7]     -5.32
##    m0[8]     -1.51
##    m0[9]     -4.96
##    m0[10]    -0.48
##    m0[11]    -5.82
##    m0[12]    -0.52
##    m0[13]    -5.40
##    m0[14]     1.70
##    m0[15]     1.38
##    m0[16]     0.81
##    m0[17]    -1.98
##    m0[18]    -6.40
##    m0[19]    -4.15
##    m0[20]     1.49
##    y0[1]      7.86
##    y0[2]      0.09
##    y0[3]      0.01
##    y0[4]      4.56
##    y0[5]      9.98
##    y0[6]      0.03
##    y0[7]      0.01
##    y0[8]      0.17
##    y0[9]      0.01
##    y0[10]     1.05
##    y0[11]     0.00
##    y0[12]     1.02
##    y0[13]     0.00
##    y0[14]     6.44
##    y0[15]     6.24
##    y0[16]     3.84
##    y0[17]     0.21
##    y0[18]     0.00
##    y0[19]     0.01
##    y0[20]     4.81
##    m1[1]      2.02
##    m1[2]      1.08
##    m1[3]      0.15
##    m1[4]     -0.79
##    m1[5]     -1.72
##    m1[6]     -2.66
##    m1[7]     -3.59
##    m1[8]     -4.53
##    m1[9]     -5.46
##    m1[10]    -6.40
##    y1[1]      7.11
##    y1[2]      3.17
##    y1[3]      0.59
##    y1[4]      1.15
##    y1[5]      0.15
##    y1[6]      0.08
##    y1[7]      0.02
##    y1[8]      0.01
##    y1[9]      0.00
##    y1[10]     0.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -8.53  -8.16 1.38 1.05 -10.95 -7.09 1.01      433      545
##    b0      9.96   9.75 1.86 1.75   7.22 13.27 1.01      597      446
##    b1     -1.98  -1.98 0.07 0.07  -2.10 -1.87 1.01      748      720
##    s       0.45   0.44 0.08 0.07   0.34  0.60 1.01      698      609
##    m0[1]   1.75   1.74 0.17 0.16   1.47  2.03 1.01      596      425
##    m0[2]  -1.68  -1.68 0.10 0.10  -1.85 -1.51 1.00     1078     1363
##    m0[3]  -4.29  -4.29 0.13 0.13  -4.50 -4.07 1.01     2155     1193
##    m0[4]   2.04   2.04 0.18 0.17   1.75  2.34 1.01      597      437
##    m0[5]   1.91   1.91 0.18 0.17   1.63  2.20 1.01      597      431
##    m0[6]  -2.68  -2.68 0.11 0.10  -2.86 -2.51 1.00     1777     1166
##    m0[7]  -5.33  -5.33 0.16 0.15  -5.59 -5.07 1.00     1917     1213
##    m0[8]  -1.50  -1.50 0.11 0.10  -1.67 -1.33 1.00      990     1303
##    m0[9]  -4.96  -4.97 0.15 0.14  -5.21 -4.72 1.01     2010     1159
##    m0[10] -0.47  -0.46 0.12 0.11  -0.66 -0.28 1.01      709      747
##    m0[11] -5.83  -5.84 0.17 0.17  -6.12 -5.55 1.00     1764     1238
##    m0[12] -0.50  -0.50 0.12 0.11  -0.70 -0.32 1.01      714      747
##    m0[13] -5.41  -5.41 0.16 0.16  -5.67 -5.15 1.00     1897     1197
##    m0[14]  1.72   1.72 0.17 0.16   1.45  2.01 1.01      597      425
##    m0[15]  1.40   1.40 0.16 0.16   1.14  1.68 1.01      601      460
##    m0[16]  0.83   0.83 0.15 0.14   0.59  1.07 1.01      613      513
##    m0[17] -1.97  -1.97 0.10 0.10  -2.15 -1.81 1.00     1259     1347
##    m0[18] -6.41  -6.42 0.19 0.19  -6.72 -6.10 1.00     1593     1167
##    m0[19] -4.16  -4.16 0.13 0.13  -4.36 -3.94 1.01     2172     1109
##    m0[20]  1.52   1.51 0.17 0.16   1.25  1.79 1.01      599      459
##    y0[1]   6.54   5.90 3.51 2.67   2.58 12.76 1.00     1532     1513
##    y0[2]   0.21   0.19 0.10 0.08   0.09  0.40 1.00     1993     1645
##    y0[3]   0.02   0.01 0.01 0.01   0.01  0.03 1.00     2013     2023
##    y0[4]   8.72   7.70 4.78 3.58   3.43 16.92 1.00     1931     1727
##    y0[5]   7.50   6.68 3.88 3.15   3.01 14.60 1.00     1528     1688
##    y0[6]   0.08   0.07 0.04 0.03   0.03  0.15 1.00     2129     1840
##    y0[7]   0.01   0.00 0.00 0.00   0.00  0.01 1.00     1877     1845
##    y0[8]   0.25   0.22 0.12 0.10   0.10  0.46 1.00     1843     1739
##    y0[9]   0.01   0.01 0.00 0.00   0.00  0.01 1.00     1684     1741
##    y0[10]  0.71   0.63 0.38 0.28   0.29  1.39 1.00     2158     1994
##    y0[11]  0.00   0.00 0.00 0.00   0.00  0.01 1.00     1939     1970
##    y0[12]  0.66   0.59 0.33 0.26   0.28  1.27 1.00     1824     2015
##    y0[13]  0.01   0.00 0.00 0.00   0.00  0.01 1.00     2110     2042
##    y0[14]  6.35   5.64 3.30 2.57   2.63 12.67 1.00     1662     1690
##    y0[15]  4.67   4.07 2.48 1.85   1.90  9.33 1.00     1438     1332
##    y0[16]  2.52   2.25 1.28 1.01   1.01  4.99 1.00     1642     1443
##    y0[17]  0.16   0.14 0.08 0.06   0.06  0.31 1.00     1983     1940
##    y0[18]  0.00   0.00 0.00 0.00   0.00  0.00 1.00     1939     1858
##    y0[19]  0.02   0.02 0.01 0.01   0.01  0.04 1.00     2195     1922
##    y0[20]  5.06   4.42 2.63 2.01   2.04  9.72 1.00     1545     1584
##    m1[1]   2.04   2.04 0.18 0.17   1.75  2.34 1.01      597      437
##    m1[2]   1.10   1.10 0.15 0.15   0.86  1.36 1.01      606      498
##    m1[3]   0.16   0.16 0.13 0.13  -0.05  0.37 1.01      647      584
##    m1[4]  -0.78  -0.77 0.11 0.11  -0.96 -0.60 1.00      762      870
##    m1[5]  -1.72  -1.72 0.10 0.10  -1.89 -1.55 1.00     1097     1379
##    m1[6]  -2.66  -2.65 0.11 0.10  -2.83 -2.48 1.00     1761     1166
##    m1[7]  -3.60  -3.59 0.12 0.12  -3.79 -3.41 1.01     2158     1211
##    m1[8]  -4.54  -4.54 0.14 0.14  -4.76 -4.31 1.01     2113     1193
##    m1[9]  -5.47  -5.48 0.16 0.16  -5.74 -5.21 1.00     1877     1197
##    m1[10] -6.41  -6.42 0.19 0.19  -6.72 -6.10 1.00     1593     1167
##    y1[1]   8.74   7.89 4.28 3.53   3.54 16.98 1.00     1838     1896
##    y1[2]   3.39   3.06 1.74 1.31   1.38  6.38 1.00     1836     1512
##    y1[3]   1.34   1.21 0.69 0.54   0.54  2.54 1.00     1546     1816
##    y1[4]   0.50   0.45 0.26 0.19   0.22  0.95 1.00     1972     1884
##    y1[5]   0.21   0.18 0.11 0.08   0.09  0.41 1.00     1827     1920
##    y1[6]   0.08   0.07 0.04 0.03   0.03  0.15 1.00     1931     1804
##    y1[7]   0.03   0.03 0.02 0.01   0.01  0.06 1.00     1983     2014
##    y1[8]   0.01   0.01 0.01 0.00   0.01  0.02 1.00     2046     1753
##    y1[9]   0.00   0.00 0.00 0.00   0.00  0.01 1.00     1707     1709
##    y1[10]  0.00   0.00 0.00 0.00   0.00  0.00 1.00     2272     1860

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,log(y),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

ex8-4

growth curve  

y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)  
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)
data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=10> b1;
  real<lower=0> s;
}
model{
  y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(1-exp(-b1*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(1-exp(-b1*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)
## Initial log joint probability = -1866.79 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       26      -9.36325   0.000146461    0.00351921      0.9112      0.9112       60    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -9.36
##    b0        11.49
##    b1         0.36
##    s          0.97
##    m0[1]      8.15
##    m0[2]      3.43
##    m0[3]      9.88
##    m0[4]     10.42
##    m0[5]      1.17
##    m0[6]     10.72
##    m0[7]      9.47
##    m0[8]      2.22
##    m0[9]      5.98
##    m0[10]    11.02
##    m0[11]     8.43
##    m0[12]    10.27
##    m0[13]    10.03
##    m0[14]     8.28
##    m0[15]     2.28
##    m0[16]    10.11
##    m0[17]    10.98
##    m0[18]     1.17
##    m0[19]     9.36
##    m0[20]    11.00
##    y0[1]      9.89
##    y0[2]      1.47
##    y0[3]      9.72
##    y0[4]     10.42
##    y0[5]      0.81
##    y0[6]     11.13
##    y0[7]      9.08
##    y0[8]      1.46
##    y0[9]      6.07
##    y0[10]    11.64
##    y0[11]     8.58
##    y0[12]    10.48
##    y0[13]    10.74
##    y0[14]     7.42
##    y0[15]     2.36
##    y0[16]     7.12
##    y0[17]    10.54
##    y0[18]     0.84
##    y0[19]    10.56
##    y0[20]    10.19
##    m1[1]      1.17
##    m1[2]      4.17
##    m1[3]      6.30
##    m1[4]      7.81
##    m1[5]      8.88
##    m1[6]      9.64
##    m1[7]     10.18
##    m1[8]     10.56
##    m1[9]     10.83
##    m1[10]    11.02
##    y1[1]      2.04
##    y1[2]      4.21
##    y1[3]      6.93
##    y1[4]      6.87
##    y1[5]      6.87
##    y1[6]      9.52
##    y1[7]     11.19
##    y1[8]      9.10
##    y1[9]     11.35
##    y1[10]    10.78
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -9.76  -9.40 1.36 1.11 -12.55 -8.29 1.00      567      552
##    b0     11.51  11.43 0.81 0.72  10.36 12.95 1.00      678      677
##    b1      0.38   0.37 0.08 0.07   0.27  0.53 1.00      598      617
##    s       1.10   1.07 0.20 0.19   0.82  1.47 1.00      871      747
##    m0[1]   8.18   8.19 0.42 0.41   7.49  8.84 1.00      890     1162
##    m0[2]   3.50   3.47 0.41 0.37   2.89  4.25 1.00      655      537
##    m0[3]   9.86   9.87 0.30 0.28   9.38 10.32 1.00     2104     1269
##    m0[4]  10.38  10.39 0.33 0.32   9.84 10.93 1.00     1382     1083
##    m0[5]   1.21   1.19 0.18 0.16   0.95  1.53 1.00      639      494
##    m0[6]  10.68  10.69 0.39 0.37  10.03 11.32 1.00      992      889
##    m0[7]   9.47   9.47 0.31 0.29   8.93  9.95 1.00     1822     1286
##    m0[8]   2.28   2.25 0.30 0.27   1.84  2.84 1.00      646      494
##    m0[9]   6.05   6.03 0.51 0.47   5.25  6.93 1.00      704      660
##    m0[10] 10.98  10.99 0.49 0.46  10.19 11.79 1.00      810      764
##    m0[11]  8.45   8.46 0.40 0.38   7.79  9.07 1.00      962     1151
##    m0[12] 10.24  10.24 0.31 0.30   9.73 10.75 1.00     1671     1028
##    m0[13] 10.00  10.01 0.30 0.28   9.52 10.48 1.00     2058     1102
##    m0[14]  8.30   8.31 0.41 0.39   7.63  8.95 1.00      920     1191
##    m0[15]  2.34   2.31 0.31 0.28   1.89  2.91 1.00      647      494
##    m0[16] 10.09  10.09 0.30 0.28   9.59 10.57 1.00     2000     1042
##    m0[17] 10.94  10.95 0.47 0.45  10.17 11.73 1.00      828      847
##    m0[18]  1.21   1.19 0.17 0.16   0.95  1.53 1.00      639      494
##    m0[19]  9.36   9.36 0.32 0.30   8.80  9.84 1.00     1650     1323
##    m0[20] 10.96  10.97 0.48 0.46  10.18 11.76 1.00      818      764
##    y0[1]   8.20   8.19 1.21 1.18   6.24 10.17 1.00     1935     1968
##    y0[2]   3.55   3.53 1.23 1.18   1.67  5.69 1.00     1642     1913
##    y0[3]   9.86   9.86 1.17 1.10   7.99 11.78 1.00     1910     1833
##    y0[4]  10.36  10.35 1.14 1.08   8.51 12.25 1.00     1937     1881
##    y0[5]   1.20   1.21 1.12 1.06  -0.67  3.01 1.00     1930     1915
##    y0[6]  10.68  10.71 1.22 1.12   8.67 12.58 1.00     1899     1786
##    y0[7]   9.47   9.48 1.14 1.06   7.63 11.30 1.00     1617     1805
##    y0[8]   2.29   2.30 1.18 1.11   0.38  4.23 1.00     1799     1394
##    y0[9]   6.07   6.05 1.24 1.19   4.08  8.17 1.00     1583     1750
##    y0[10] 10.98  11.01 1.25 1.17   8.94 13.05 1.00     1504     1824
##    y0[11]  8.43   8.45 1.17 1.17   6.55 10.32 1.00     1727     1958
##    y0[12] 10.22  10.22 1.16 1.10   8.27 12.09 1.00     1785     1641
##    y0[13] 10.01  10.00 1.16 1.12   8.12 11.99 1.00     2165     1954
##    y0[14]  8.29   8.28 1.16 1.09   6.37 10.17 1.00     1838     1856
##    y0[15]  2.33   2.28 1.15 1.14   0.52  4.33 1.00     1686     1891
##    y0[16] 10.07  10.07 1.16 1.13   8.11 11.94 1.00     1910     1956
##    y0[17] 10.95  11.00 1.23 1.21   8.90 12.92 1.00     1938     1894
##    y0[18]  1.21   1.20 1.16 1.08  -0.71  3.07 1.00     1886     1820
##    y0[19]  9.32   9.32 1.17 1.10   7.38 11.30 1.00     1980     1933
##    y0[20] 10.94  10.93 1.23 1.21   8.96 12.94 1.00     1611     1514
##    m1[1]   1.21   1.19 0.17 0.16   0.95  1.53 1.00      639      494
##    m1[2]   4.25   4.22 0.46 0.42   3.56  5.07 1.00      661      527
##    m1[3]   6.37   6.35 0.50 0.47   5.57  7.23 1.00      714      725
##    m1[4]   7.85   7.84 0.44 0.43   7.12  8.56 1.00      826      951
##    m1[5]   8.89   8.90 0.36 0.33   8.27  9.45 1.00     1162     1341
##    m1[6]   9.62   9.63 0.30 0.28   9.10 10.10 1.00     1993     1323
##    m1[7]  10.15  10.15 0.31 0.29   9.65 10.64 1.00     1865     1164
##    m1[8]  10.52  10.52 0.36 0.34   9.94 11.09 1.00     1171      928
##    m1[9]  10.79  10.80 0.42 0.40  10.11 11.48 1.00      908      859
##    m1[10] 10.98  10.99 0.49 0.46  10.19 11.79 1.00      810      764
##    y1[1]   1.23   1.22 1.14 1.12  -0.62  3.13 1.00     1877     1842
##    y1[2]   4.25   4.27 1.19 1.11   2.34  6.21 1.00     1748     1522
##    y1[3]   6.36   6.35 1.26 1.19   4.34  8.41 1.00     1601     1772
##    y1[4]   7.86   7.85 1.22 1.17   5.87  9.82 1.00     1853     1575
##    y1[5]   8.89   8.90 1.14 1.06   7.00 10.68 1.00     2045     1900
##    y1[6]   9.63   9.66 1.17 1.16   7.77 11.49 1.00     1930     1772
##    y1[7]  10.18  10.18 1.16 1.09   8.33 12.02 1.00     2176     1897
##    y1[8]  10.51  10.54 1.19 1.09   8.55 12.42 1.00     1801     1856
##    y1[9]  10.80  10.79 1.16 1.12   8.93 12.71 1.00     1849      991
##    y1[10] 11.01  11.00 1.19 1.15   9.08 12.99 1.00     1692     2012

ex8-5

sigmoid curve

y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)
data{
  int N;
  vector[N] x;
  int Ym;
  array[N] int y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=100> b1;
}
model{
  y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
  array[N] int y0;
  for(i in 1:N){
    y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
  }
  array[N1] int y1;
  for(i in 1:N1){
    y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
  }
}
n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-5.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "y0"   "y1"
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.61 -16.31 1.03 0.72 -18.66 -15.65 1.01      626      584
##    b0       4.00   3.98 0.26 0.26   3.59   4.44 1.00     1191     1273
##    b1       2.08   2.01 0.47 0.43   1.45   2.92 1.01      378      257
##    y0[1]    0.09   0.00 0.31 0.00   0.00   1.00 1.00     2023     2027
##    y0[2]    2.68   3.00 1.61 1.48   0.00   6.00 1.00     1884     2041
##    y0[3]   10.00  10.00 0.06 0.00  10.00  10.00 1.00     2030       NA
##    y0[4]   10.00  10.00 0.07 0.00  10.00  10.00 1.00     1700       NA
##    y0[5]    9.97  10.00 0.16 0.00  10.00  10.00 1.00     1681       NA
##    y0[6]    0.31   0.00 0.61 0.00   0.00   2.00 1.00     1500     1462
##    y0[7]    0.25   0.00 0.53 0.00   0.00   1.00 1.00     1797     1824
##    y0[8]    9.93  10.00 0.28 0.00   9.00  10.00 1.00     1589       NA
##    y0[9]    0.24   0.00 0.50 0.00   0.00   1.00 1.00     1707     1754
##    y0[10]   9.83  10.00 0.43 0.00   9.00  10.00 1.00     1491       NA
##    y0[11]  10.00  10.00 0.06 0.00  10.00  10.00 1.00     2026       NA
##    y0[12]   9.96  10.00 0.20 0.00  10.00  10.00 1.00     1607       NA
##    y0[13]   9.96  10.00 0.20 0.00  10.00  10.00 1.00     1922       NA
##    y0[14]   9.84  10.00 0.44 0.00   9.00  10.00 1.00      955       NA
##    y0[15]   1.87   2.00 1.41 1.48   0.00   4.00 1.00     1913     1915
##    y0[16]  10.00  10.00 0.04 0.00  10.00  10.00 1.00     2023       NA
##    y0[17]   0.07   0.00 0.27 0.00   0.00   1.00 1.00     1576     1599
##    y0[18]   9.90  10.00 0.33 0.00   9.00  10.00 1.00     1440       NA
##    y0[19]  10.00  10.00 0.06 0.00  10.00  10.00 1.00     1655       NA
##    y0[20]   0.01   0.00 0.11 0.00   0.00   0.00 1.00     1744     1744
##    y1[1]    0.01   0.00 0.08 0.00   0.00   0.00 1.00     2043     2043
##    y1[2]    0.05   0.00 0.22 0.00   0.00   0.00 1.00     1747     1777
##    y1[3]    0.26   0.00 0.53 0.00   0.00   1.00 1.00     1739     1711
##    y1[4]    1.43   1.00 1.21 1.48   0.00   4.00 1.00     1976     1785
##    y1[5]    5.18   5.00 1.94 1.48   2.00   8.00 1.00     1548     1548
##    y1[6]    8.59   9.00 1.38 1.48   6.00  10.00 1.00     1011       NA
##    y1[7]    9.69  10.00 0.63 0.00   8.00  10.00 1.00     1257       NA
##    y1[8]    9.93  10.00 0.27 0.00   9.00  10.00 1.00     1600       NA
##    y1[9]    9.98  10.00 0.13 0.00  10.00  10.00 1.00     1861       NA
##    y1[10]  10.00  10.00 0.06 0.00  10.00  10.00 1.00     2029       NA

y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=ymed),xy,col='black')

ex8-6

up and down

y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))
data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=200> b0;
  real<lower=0,upper=1> b1;
  real<lower=0,upper=1> b2;
  real<lower=0,upper=10> s;
}
model{
  y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-6.stan')

fn(mdl,data)
## Initial log joint probability = -1235.08 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       69       3.80917   1.97211e-05     0.0160704           1           1      117    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__       3.81
##    b0       104.77
##    b1         0.03
##    b2         0.19
##    s          0.50
##    m0[1]     57.76
##    m0[2]     60.85
##    m0[3]     38.52
##    m0[4]     60.80
##    m0[5]     32.73
##    m0[6]     46.81
##    m0[7]     61.40
##    m0[8]     61.41
##    m0[9]     22.94
##    m0[10]    28.87
##    m0[11]    37.08
##    m0[12]    39.17
##    m0[13]    24.72
##    m0[14]    49.48
##    m0[15]    23.99
##    m0[16]    57.04
##    m0[17]    54.01
##    m0[18]    21.98
##    m0[19]    44.47
##    m0[20]    60.73
##    y0[1]     57.69
##    y0[2]     60.67
##    y0[3]     38.61
##    y0[4]     61.40
##    y0[5]     32.52
##    y0[6]     46.53
##    y0[7]     61.76
##    y0[8]     61.27
##    y0[9]     22.87
##    y0[10]    28.73
##    y0[11]    37.31
##    y0[12]    38.96
##    y0[13]    23.74
##    y0[14]    49.43
##    y0[15]    24.21
##    y0[16]    56.89
##    y0[17]    54.50
##    y0[18]    21.89
##    y0[19]    44.42
##    y0[20]    60.66
##    m1[1]     28.87
##    m1[2]     57.76
##    m1[3]     61.06
##    m1[4]     56.14
##    m1[5]     49.17
##    m1[6]     42.24
##    m1[7]     36.00
##    m1[8]     30.57
##    m1[9]     25.93
##    m1[10]    21.98
##    y1[1]     29.51
##    y1[2]     57.63
##    y1[3]     61.11
##    y1[4]     56.19
##    y1[5]     48.70
##    y1[6]     42.11
##    y1[7]     36.03
##    y1[8]     31.22
##    y1[9]     26.03
##    y1[10]    21.45
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__    -0.61  -0.24 1.66 1.44  -3.63   1.27 1.01      363      311
##    b0     104.90 104.81 1.81 1.73 102.10 108.02 1.01     1072      800
##    b1       0.03   0.03 0.00 0.00   0.03   0.03 1.01     1148      972
##    b2       0.19   0.19 0.00 0.00   0.19   0.20 1.01     1056      868
##    s        0.60   0.58 0.12 0.10   0.46   0.80 1.03      203      199
##    m0[1]   57.76  57.76 0.22 0.21  57.41  58.12 1.00     1434     1357
##    m0[2]   60.86  60.86 0.21 0.21  60.53  61.19 1.00     1762     1443
##    m0[3]   38.52  38.50 0.37 0.35  37.94  39.12 1.01     1226     1296
##    m0[4]   60.80  60.80 0.21 0.21  60.47  61.13 1.00     1745     1443
##    m0[5]   32.72  32.72 0.20 0.19  32.40  33.05 1.01     1845     1305
##    m0[6]   46.81  46.81 0.19 0.18  46.49  47.13 1.00     1623     1217
##    m0[7]   61.41  61.40 0.21 0.21  61.06  61.75 1.00     1959     1478
##    m0[8]   61.41  61.41 0.21 0.21  61.07  61.75 1.00     1952     1490
##    m0[9]   22.93  22.93 0.25 0.24  22.52  23.33 1.01     1457     1228
##    m0[10]  28.87  28.85 0.32 0.31  28.36  29.39 1.01     1185     1221
##    m0[11]  37.07  37.07 0.19 0.17  36.76  37.38 1.00     1932     1245
##    m0[12]  39.17  39.17 0.18 0.17  38.86  39.48 1.00     1911     1251
##    m0[13]  24.71  24.71 0.24 0.24  24.31  25.10 1.01     1512     1275
##    m0[14]  49.48  49.47 0.20 0.19  49.15  49.82 1.00     1509     1161
##    m0[15]  23.98  23.98 0.25 0.24  23.58  24.37 1.01     1488     1279
##    m0[16]  57.04  57.03 0.30 0.28  56.56  57.54 1.00     1634     1237
##    m0[17]  54.01  54.01 0.22 0.21  53.66  54.37 1.00     1397     1266
##    m0[18]  21.96  21.97 0.26 0.24  21.55  22.37 1.01     1431     1227
##    m0[19]  44.47  44.47 0.19 0.17  44.15  44.77 1.00     1733     1317
##    m0[20]  60.73  60.73 0.21 0.21  60.39  61.06 1.00     1731     1431
##    y0[1]   57.77  57.79 0.65 0.61  56.71  58.84 1.00     1883     1838
##    y0[2]   60.86  60.84 0.62 0.58  59.86  61.87 1.00     1791     1629
##    y0[3]   38.53  38.51 0.72 0.67  37.39  39.69 1.01     1756     1540
##    y0[4]   60.80  60.80 0.63 0.60  59.76  61.80 1.00     1816     1810
##    y0[5]   32.75  32.73 0.64 0.62  31.70  33.83 1.00     2045     1774
##    y0[6]   46.81  46.80 0.65 0.62  45.73  47.83 1.00     1818     1847
##    y0[7]   61.41  61.40 0.64 0.61  60.37  62.51 1.00     2014     1434
##    y0[8]   61.38  61.37 0.65 0.61  60.34  62.45 1.00     2002     1780
##    y0[9]   22.91  22.90 0.68 0.61  21.83  24.02 1.00     1894     1497
##    y0[10]  28.86  28.85 0.71 0.67  27.69  30.01 1.00     1545     1772
##    y0[11]  37.06  37.06 0.64 0.62  36.00  38.05 1.00     2076     1829
##    y0[12]  39.16  39.16 0.65 0.62  38.11  40.19 1.00     1858     1756
##    y0[13]  24.68  24.69 0.67 0.62  23.56  25.78 1.00     2124     1587
##    y0[14]  49.47  49.45 0.64 0.61  48.45  50.52 1.00     1831     1972
##    y0[15]  23.99  23.98 0.66 0.63  22.94  25.07 1.00     1698     1654
##    y0[16]  57.04  57.04 0.69 0.66  55.94  58.16 1.00     1964     1970
##    y0[17]  54.02  54.01 0.64 0.64  52.99  55.10 1.00     1978     1902
##    y0[18]  21.97  21.97 0.69 0.66  20.85  23.06 1.00     1853     1700
##    y0[19]  44.47  44.47 0.65 0.62  43.41  45.54 1.00     2133     1659
##    y0[20]  60.72  60.70 0.66 0.64  59.66  61.80 1.00     1942     1701
##    m1[1]   28.87  28.85 0.32 0.31  28.36  29.39 1.01     1185     1221
##    m1[2]   57.76  57.74 0.29 0.27  57.29  58.23 1.00     1675     1279
##    m1[3]   61.06  61.06 0.21 0.20  60.73  61.39 1.00     1802     1452
##    m1[4]   56.15  56.14 0.22 0.21  55.79  56.51 1.00     1400     1299
##    m1[5]   49.17  49.17 0.20 0.19  48.84  49.50 1.00     1521     1161
##    m1[6]   42.24  42.23 0.18 0.17  41.93  42.54 1.00     1821     1312
##    m1[7]   35.99  35.99 0.19 0.18  35.68  36.30 1.00     1923     1187
##    m1[8]   30.56  30.56 0.21 0.21  30.23  30.91 1.01     1772     1362
##    m1[9]   25.92  25.92 0.24 0.23  25.54  26.30 1.01     1559     1279
##    m1[10]  21.96  21.97 0.26 0.24  21.55  22.37 1.01     1431     1227
##    y1[1]   28.86  28.85 0.71 0.68  27.74  29.99 1.00     1794     1529
##    y1[2]   57.77  57.77 0.66 0.64  56.71  58.84 1.00     1984     1519
##    y1[3]   61.07  61.05 0.64 0.61  59.98  62.11 1.00     2060     1563
##    y1[4]   56.14  56.13 0.66 0.62  55.07  57.21 1.00     1919     1341
##    y1[5]   49.18  49.21 0.65 0.62  48.09  50.20 1.00     2206     1883
##    y1[6]   42.23  42.23 0.64 0.59  41.16  43.27 1.00     1941     1507
##    y1[7]   36.03  36.02 0.63 0.61  35.00  37.06 1.00     1995     1914
##    y1[8]   30.57  30.56 0.65 0.62  29.52  31.58 1.01     2077     1832
##    y1[9]   25.92  25.94 0.65 0.65  24.86  26.97 1.00     1988     1584
##    y1[10]  21.97  21.97 0.66 0.63  20.88  23.07 1.00     2006     1938